Parámetros de cuenca con r.basin

library(rgrass7)
## Loading required package: XML
## GRASS GIS interface loaded with GRASS version: GRASS 7.8.2 (2019)
## and location: rdom
gisdbase <- 'grass-data-test' #Base de datos de GRASS GIS
wd <- getwd() #Directorio de trabajo
wd
## [1] "/home/edelte/unidad-0-asignacion-99-mi-manuscrito-edeltejeda"
loc <- initGRASS(gisBase = "/usr/lib/grass78/",
                 home = wd,
                 gisDbase = paste(wd, gisdbase, sep = '/'),
                 location = 'c_ozama',
                 mapset = "PERMANENT",
                 override = TRUE)
gmeta()
## gisdbase    /home/edelte/unidad-0-asignacion-99-mi-manuscrito-edeltejeda/grass-data-test 
## location    c_ozama 
## mapset      PERMANENT 
## rows        1 
## columns     1 
## north       1 
## south       0 
## west        0 
## east        1 
## nsres       1 
## ewres       1 
## projection  NA
execGRASS(
  'g.list',
  flags = 't',
  parameters = list(
    type = c('raster', 'vector')
  )
)

Convertir a números enteros la extensión y la resolución del DEM

library(raster)
## Loading required package: sp
rutadem <- 'data/dem.tif'
rawextent <- extent(raster(rutadem))
rawextent
## class      : Extent 
## xmin       : 314490.3 
## xmax       : 460848.8 
## ymin       : 2008160 
## ymax       : 2103177
devtools::source_url('https://raw.githubusercontent.com/geofis/rgrass/master/integerextent.R')
## SHA-1 hash of file is 65b57fe7cb20cd1976572cd0a7e98def9e95d83c
devtools::source_url('https://raw.githubusercontent.com/geofis/rgrass/master/xyvector.R')
## SHA-1 hash of file is 89fa5ae436d9a7d7a0c799b789c560eb5e421cfd
newextent <- intext(e = rawextent, r = 90, type = 'inner')
newextent
## class      : Extent 
## xmin       : 314550 
## xmax       : 460800 
## ymin       : 2008170 
## ymax       : 2103120
gdalUtils::gdalwarp(
  srcfile = 'data/dem.tif',
  dstfile = 'data/demint.tif',
  te = xyvector(newextent),
  tr = c(90,90),
  r = 'bilinear',
  overwrite = T
)
## Registered S3 method overwritten by 'R.oo':
##   method        from       
##   throw.default R.methodsS3
## NULL

Importar a sesión de GRASS

rutademint <- 'data/demint.tif'
execGRASS(
  "g.proj",
  flags = c('t','c'),
  georef=rutademint)
gmeta()
execGRASS(
  "r.in.gdal",
  flags='overwrite',
  parameters=list(
    input=rutademint,
    output="demint"
  )
)
execGRASS(
  "g.region",
  parameters=list(
    raster = "demint",
    align = "demint"
  )
)
gmeta()
## gisdbase    /home/edelte/unidad-0-asignacion-99-mi-manuscrito-edeltejeda/grass-data-test 
## location    c_ozama 
## mapset      PERMANENT 
## rows        1055 
## columns     1625 
## north       2103120 
## south       2008170 
## west        314550 
## east        460800 
## nsres       90 
## ewres       90 
## projection  +proj=utm +no_defs +zone=19 +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000 +type=crs +to_meter=1
execGRASS(
  'g.list',
  flags = 't',
  parameters = list(
    type = c('raster', 'vector')
  )
)
## raster/demint

Generar red de drenaje para obtener coordenada posteriormente

execGRASS(
  "r.stream.extract",
  flags = c('overwrite','quiet'),
  parameters = list(
    elevation = 'demint',
    threshold = 80,
    stream_raster = 'stream-de-rstr',
    stream_vector = 'stream_de_rstr'
  )
)
execGRASS(
  'g.list',
  flags = 't',
  parameters = list(
    type = c('raster', 'vector')
  )
)
## raster/demint
## raster/stream-de-rstr
## vector/stream_de_rstr

Obtener coordenada

library(sp)
use_sp()
library(mapview)
netw <- spTransform(
  readVECT('stream_de_rstr'),
  CRSobj = CRS("+init=epsg:4326"))
mapview(netw, col.regions = 'blue', legend = FALSE)
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'

Transformar coordenada a EPSG:32619 como número entero

source('my-trans.R')
outlet <- as.integer(my_trans(c(-69.88117,18.47503)))

Ejecutar r.basin

pref <- 'rbasin_ozama'
execGRASS(
  "r.basin",
  flags = 'overwrite',
  parameters = list(
    map = 'demint',
    prefix = pref,
    coordinates = outlet,
    threshold = 80,
    dir = 'salidas-rbasin/ozama'
  )
)
execGRASS(
  'g.list',
  flags = 't',
  parameters = list(
    type = c('raster', 'vector')
  )
)
## raster/demint
## raster/rbasin_ozama_demint_accumulation
## raster/rbasin_ozama_demint_aspect
## raster/rbasin_ozama_demint_dist2out
## raster/rbasin_ozama_demint_drainage
## raster/rbasin_ozama_demint_hack
## raster/rbasin_ozama_demint_hillslope_distance
## raster/rbasin_ozama_demint_horton
## raster/rbasin_ozama_demint_shreve
## raster/rbasin_ozama_demint_slope
## raster/rbasin_ozama_demint_strahler
## raster/stream-de-rstr
## vector/rbasin_ozama_demint_basin
## vector/rbasin_ozama_demint_mainchannel
## vector/rbasin_ozama_demint_network
## vector/rbasin_ozama_demint_outlet
## vector/rbasin_ozama_demint_outlet_snap
## vector/stream_de_rstr

Si r.basin arrojara error (sólo en el caso de error, no en caso de advertencia), ejecutar este bloque para borrar las salidas anteriores y reejecutar el r.basin:

execGRASS(
  "g.remove",
  flags = 'f',
  parameters = list(
    type = c('raster','vector'),
    pattern = paste0(pref, '*')
  )
)

Cargar los vectoriales transformados a EPSG:4326 para visualizar en leaflet

rbnetw <- spTransform(
  readVECT('rbasin_ozama_demint_network'),
  CRSobj = CRS("+init=epsg:4326"))
rbnetw
rbmain <- spTransform(
  readVECT('rbasin_ozama_demint_mainchannel'),
  CRSobj = CRS("+init=epsg:4326"))
rbmain
rbbasin <- spTransform(
  readVECT('rbasin_ozama_demint_basin'),
  CRSobj = CRS("+init=epsg:4326"))
rbbasin
library(leaflet)
leaflet() %>%
  addProviderTiles(providers$Stamen.Terrain, group = 'terrain') %>%
  addPolylines(data = rbnetw, weight = 3, opacity = 0.7) %>% 
  addPolylines(data = rbmain, weight = 3, opacity = 0.7, color = 'red') %>% 
  addPolygons(data = rbbasin) %>% 
  leafem::addHomeButton(extent(rbbasin), 'Ver cuenca')

Explorar los parámetros de cuenca

library(readr)
rbozamapar1 <- read_csv("salidas-rbasin/ozama/rbasin_ozama_demint_parametersT.csv")
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   Rectangle_containing_basin_N_W = col_character(),
##   Rectangle_containing_basin_S_E = col_character()
## )
## See spec(...) for full column specifications.
rbozamapar1 %>% tibble::as_tibble()
## # A tibble: 1 x 34
##        x      y Easting_Centroi… Northing_Centro… Rectangle_conta…
##    <dbl>  <dbl>            <dbl>            <dbl> <chr>           
## 1 406935 2.04e6           413955          2069865 ('367560', '209…
## # … with 29 more variables: Rectangle_containing_basin_S_E <chr>,
## #   Area_of_basin_km2 <dbl>, Perimeter_of_basin_km <dbl>,
## #   Max_Elevation <dbl>, Min_Elevation <dbl>, Elevation_Difference <dbl>,
## #   Mean_Elevation <dbl>, Mean_Slope <dbl>,
## #   Length_of_Directing_Vector_km <dbl>,
## #   Prevalent_Orientation_deg_from_north_ccw <dbl>,
## #   Compactness_Coefficient <dbl>, Circularity_Ratio <dbl>,
## #   Topological_Diameter <dbl>, Elongation_Ratio <dbl>,
## #   Shape_Factor <dbl>, Concentration_Time_hr <dbl>,
## #   Length_of_Mainchannel_km <dbl>,
## #   Mean_slope_of_mainchannel_percent <dbl>,
## #   Mean_hillslope_length_m <dbl>, Magnitudo <dbl>,
## #   Max_order_Strahler <dbl>, Number_of_streams <dbl>,
## #   Total_Stream_Length_km <dbl>, First_order_stream_frequency <dbl>,
## #   Drainage_Density_km_over_km2 <dbl>, Bifurcation_Ratio_Horton <dbl>,
## #   Length_Ratio_Horton <dbl>, Area_ratio_Horton <dbl>,
## #   Slope_ratio_Horton <dbl>
rbozamapar2 <- read_csv(
  "salidas-rbasin/ozama/rbasin_ozama_demint_parameters.csv",
  skip=2, col_names = c('Parameter', 'Value'))
## Parsed with column specification:
## cols(
##   Parameter = col_character(),
##   Value = col_character()
## )
rbozamapar2 %>% print(n=Inf)
## # A tibble: 32 x 2
##    Parameter                                             Value             
##    <chr>                                                 <chr>             
##  1 Easting Centroid of basin                             413955.00         
##  2 Northing Centroid of basin                            2069865.00        
##  3 Rectangle containing basin N-W                        ('367560', '20983…
##  4 Rectangle containing basin S-E                        ('455220', '20399…
##  5 Area of basin [km^2]                                  3408.388875       
##  6 Perimeter of basin [km]                               358.481006282309  
##  7 Max Elevation [m s.l.m.]                              913.3936          
##  8 Min Elevation [m s.l.m.]                              -3.278881         
##  9 Elevation Difference [m]                              916.672481        
## 10 Mean Elevation                                        106.245           
## 11 Mean Slope                                            3.47              
## 12 Length of Directing Vector [km]                       27.788971913332816
## 13 Prevalent Orientation [degree from north, counterclo… 1.3166010130150796
## 14 Compactness Coefficient                               5.441724127967625 
## 15 Circularity Ratio                                     0.333293391887047…
## 16 Topological Diameter                                  168.0             
## 17 Elongation Ratio                                      0.5008280978139896
## 18 Shape Factor                                          25.91243324088714 
## 19 Concentration Time (Giandotti, 1934) [hr]             17.787167537498174
## 20 Length of Mainchannel [km]                            131.534883016     
## 21 Mean slope of mainchannel [percent]                   1.1417035909306994
## 22 Mean hillslope length [m]                             4962.319          
## 23 Magnitudo                                             1096.0            
## 24 Max order (Strahler)                                  7                 
## 25 Number of streams                                     1615              
## 26 Total Stream Length [km]                              3117.2892         
## 27 First order stream frequency                          0.321559552091895…
## 28 Drainage Density [km/km^2]                            0.9145931741723573
## 29 Bifurcation Ratio (Horton)                            3.3679            
## 30 Length Ratio (Horton)                                 1.6319            
## 31 Area ratio (Horton)                                   3.6528            
## 32 Slope ratio (Horton)                                  1.3764

Limpiar archivo de bloqueo del conjunto de mapas de GRASS

unlink_.gislock()